Patterns in mineral-associated organic matter across scales: untangling the influence of climate, mineralogy and mycorrhizal fungi on soil carbon stabilization

#Assessing potential influence of unusual values in final data set

hist(final_data$mg.C.per.g.soil_FLF) #very high concentration of FLF C in one sample

hist(final_data$mg.C.per.g.soil_OLF)

hist(final_data$mg.C.per.g.soil_HF)

hist(final_data$mg.N.per.g.soil_FLF) #very high concentration of FLF N in one sample

hist(final_data$mg.N.per.g.soil_OLF)

hist(final_data$mg.N.per.g.soil_HF)

hist(final_data$proportion.C_FLF) 

hist(final_data$proportion.C_OLF)

hist(final_data$proportion.C_HF)

hist(final_data$proportion.N_FLF) 

hist(final_data$proportion.N_OLF)

hist(final_data$proportion.N_HF)

#removing HARV_037 from final dataset


mean_FLF_conc.C=mean(final_data$mg.C.per.g.soil_FLF) #2.09
std_dev_FLF_conc.C=sd(final_data$mg.C.per.g.soil_FLF) #1.98
#HARV_037 is 47.33, greater than 20 standard deviations above mean

mean_FLF_conc.N=mean(final_data$mg.N.per.g.soil_FLF) #0.106
std_dev_FLF_conc.N=sd(final_data$mg.N.per.g.soil_FLF) #0.205
#HARV_037 has 1.42 mg FLF N per g soil, greater than 6 standard deviations above mean

final_data=final_data %>% 
  filter(plotID != "HARV_037") 

Site Characteristics

Site	MAT (C)	MAP(mm)	Dominant tree species	Soil % clay (Mean +/- SD)	Soil Fe (ox-extractable)
LENO	18.1	1386	American Sweetgum, American Hornbeam, Possumhaw	44.7 +- 13.8	1.0 +- 0.4
DELA	17.6	1372	Water oak, Red maple, Sugarberry	37.0 +- 6.3	0.88 +- 0.2
ORNL	14.4	1340	Red maple, Sour gum, Chestnut oak	22.0 +- 12.7	0.13 +- 0.05
SERC	13.6	1075	Tulip poplar, American Beech, American Sweetgum	17.4 +- 6.2	0.27 +- 0.1
HARV	7.4	1199	Eastern Hemlock, Northern Red Oak, American Beech	5.5 +- 1.9	0.38 +- 0.1
BART	6.2	1325	American Beech, Eastern Hemlock, Red maple	3.6 +- 0.7	0.76 +- 0.5
TREE	4.8	797	Sugar maple, Red maple, Gray alder	5.4 +- 1.1	0.34 +- 0.09

	proportion MAOM C						[MAOM C]						proportion MAOM N							[MAOM N]
Predictors	Estimates	std. Error	std. Beta	standardized std. Error	p	df	Estimates	std. Error	std. Beta	standardized std. Error	p	df	Estimates	std. Error	std. Beta	standardized std. Error	p	std. p	df	Estimates	std. Error	std. Beta	standardized std. Error	p	df
Intercept	0.80	0.04	0.00	0.15	<0.001	41.00	30.15	7.53	-0.00	0.24	<0.001	41.00	0.88	0.04	0.02	0.11	<0.001	0.847	40.00	1.18	0.34	-0.00	0.18	0.001	41.00
Mycorrhizal dominance	-0.10	0.03	-0.42	0.13	0.002	41.00	-1.96	3.05	-0.07	0.10	0.524	41.00	-0.19	0.07	-0.23	0.12	0.010	0.061	40.00	-0.14	0.19	-0.08	0.11	0.453	41.00
CDI	0.15	1.23	0.02	0.18	0.905	41.00	-583.93	212.12	-0.69	0.25	0.009	41.00	0.72	1.16	0.50	0.13	0.536	0.001	40.00	-15.46	9.69	-0.32	0.20	0.118	41.00
FeOx	0.01	0.00	0.50	0.16	0.003	41.00	1.36	0.33	0.62	0.15	<0.001	41.00	0.00	0.00	0.09	0.14	0.517	0.517	40.00	0.11	0.02	0.88	0.15	<0.001	41.00
Mycorrhizal dominance * CDI													4.33	1.99	0.24	0.11	0.036	0.036	40.00
Random Effects
σ2	0.00						36.91						0.00							0.14
τ00	0.00 Site						27.96 Site						0.00 Site							0.04 Site
ICC	0.09						0.43													0.24
N	7 Site						7 Site						7 Site							7 Site
Observations	47						47						47							47
Marginal R2 / Conditional R2	0.330 / 0.391						0.359 / 0.635						0.429 / NA							0.474 / 0.602

	proportion oPOM C						[oPOM C]						proportion oPOM N							[oPOM N]
Predictors	Estimates	std. Error	std. Beta	standardized std. Error	p	df	Estimates	std. Error	std. Beta	standardized std. Error	p	df	Estimates	std. Error	std. Beta	standardized std. Error	p	std. p	df	Estimates	std. Error	std. Beta	standardized std. Error	p	df
Intercept	0.06	0.03	0.00	0.18	0.039	41.00	3.55	1.33	0.00	0.19	0.011	41.00	0.03	0.03	-0.03	0.13	0.227	0.847	40.00	0.13	0.04	0.00	0.20	0.005	41.00
Mycorrhizal dominance	0.05	0.02	0.33	0.15	0.027	41.00	1.25	0.80	0.20	0.13	0.127	41.00	0.11	0.04	0.18	0.14	0.016	0.191	40.00	0.02	0.03	0.09	0.13	0.494	41.00
CDI	0.48	0.84	0.12	0.20	0.574	41.00	-79.58	37.57	-0.44	0.21	0.040	41.00	0.10	0.69	-0.40	0.16	0.881	0.014	40.00	-2.62	1.20	-0.48	0.22	0.035	41.00
FeOx	-0.00	0.00	-0.17	0.18	0.347	41.00	0.12	0.08	0.25	0.17	0.148	41.00	0.00	0.00	0.09	0.16	0.590	0.590	40.00	0.00	0.00	0.30	0.18	0.099	41.00
Mycorrhizal dominance * CDI													-2.56	1.18	-0.28	0.13	0.036	0.036	40.00
Random Effects
σ2	0.00						2.63						0.00							0.00
τ00	0.00 Site						0.56 Site						0.00 Site							0.00 Site
ICC	0.10						0.18													0.19
N	7 Site						7 Site						7 Site							7 Site
Observations	47						47						47							47
Marginal R2 / Conditional R2	0.101 / 0.190						0.215 / 0.353						0.246 / NA							0.189 / 0.346

	proportion.C_FLF						mg.C.per.g.soil_FLF						proportion.N_FLF						mg.N.per.g.soil_FLF
Predictors	Estimates	std. Error	std. Beta	standardized std. Error	p	df	Estimates	std. Error	std. Beta	standardized std. Error	p	df	Estimates	std. Error	std. Beta	standardized std. Error	p	df	Estimates	std. Error	std. Beta	standardized std. Error	p	df
(Intercept)	0.137	0.029	-0.001	0.132	<0.001	41.000	4.720	1.502	-0.009	0.216	0.003	41.000	0.119	0.021	0.000	0.122	<0.001	41.000	0.174	0.051	-0.009	0.211	0.001	41.000
E	0.051	0.023	0.280	0.127	0.033	41.000	1.011	0.745	0.162	0.119	0.182	41.000	0.027	0.018	0.198	0.128	0.129	41.000	0.019	0.027	0.091	0.124	0.471	41.000
CDI	-0.715	0.809	-0.136	0.154	0.382	41.000	-87.344	42.376	-0.480	0.233	0.046	41.000	-1.610	0.577	-0.400	0.143	0.008	41.000	-3.068	1.429	-0.491	0.229	0.038	41.000
FeOx	-0.007	0.002	-0.498	0.150	0.002	41.000	-0.021	0.078	-0.044	0.164	0.791	41.000	-0.003	0.002	-0.241	0.145	0.105	41.000	0.000	0.003	0.008	0.170	0.965	41.000
Random Effects
σ2	0.002						2.223						0.001						0.003
τ00	0.000 Site						0.951 Site						0.000 Site						0.001 Site
ICC	0.031						0.300												0.261
N	7 Site						7 Site						7 Site						7 Site
Observations	47						47						47						47
Marginal R2 / Conditional R2	0.346 / 0.366						0.265 / 0.485						0.333 / NA						0.233 / 0.433

#what % of total C and N was MAOM C in each site?
pct.MAOM.C=final_data %>% 
  group_by(Site) %>% 
  mutate(mean.C.pro.maom=mean(proportion.C_HF),
         se.pro.maom=sd(proportion.C_HF)/sqrt(46)) %>% 
  dplyr::select(plotID,mean.C.pro.maom, se.pro.maom)

## Adding missing grouping variables: `Site`

pct.MAOM.N=final_data %>% 
    group_by(Site) %>% 
  mutate(mean.N.pro.maom=mean(proportion.N_HF),
         se.pro.maom.N=sd(proportion.N_HF)/sqrt(46)) %>% 
   dplyr::select(plotID,mean.N.pro.maom, se.pro.maom.N)

## Adding missing grouping variables: `Site`

# what % of total C and N was MAOM at the lowest %ECM and at the highest %ECM in each site?

end_members_maom_prop=final_data %>% 
  select(Site,Plot.x,E, proportion.C_HF, proportion.N_HF) %>%
  arrange(Site,E) %>% 
   group_by(Site) %>% 
  mutate(myc_rank=order(E, decreasing=T)) %>% 
  mutate(max_A=max(myc_rank)) %>% 
  filter(myc_rank==max_A) %>% 
  ungroup() %>% 
  mutate(mean_prop_C=mean(proportion.C_HF),
         mean_prop_N=mean(proportion.N_HF))

 #At max ECM dominance, avg. proportion C in MAOM was 0.75831, N was 0.85565

#At max AM dominance, avg proportion C in MAOM was 0.83787, N was 0.89959

#from AM to ECM, MAOM C prop dropped 7.956%, N dropped 4.39%

concC.fe.hf=ggplot(final_data, aes(x=FeOx, y=mg.C.per.g.soil_HF, color=Site))+
   geom_smooth(method="lm",size=2, se=F)+
  geom_point()+
   scale_colour_viridis_d(option="C")+
  labs(x=" FeOx (mg/g soil)", y="[MAOM C] (mg/g soil)", colour="Site")+
  theme_cowplot()+
  theme(legend.title=element_text(hjust=0.5, size=16), legend.text=element_text(size=15), axis.text=element_text(size=15), axis.title=element_text(size=19))

propC.fe.hf=ggplot(final_data, aes(x=FeOx, y=proportion.C_HF, color=Site))+
   geom_smooth(method="lm",size=2, se=F)+
  geom_point()+
   scale_colour_viridis_d(option="C")+
  labs(x=" FeOx (mg/g soil)", y=expression("proportion C"[MAOM]), colour="Site")+
  theme_cowplot()+
  theme(legend.title=element_text(hjust=0.5, size=16), legend.text=element_text(size=15), axis.text=element_text(size=15), axis.title=element_text(size=19))


concN.fe.hf=ggplot(final_data, aes(x=FeOx, y=mg.N.per.g.soil_HF, color=Site))+
   geom_smooth(method="lm",size=2, se=F)+
  geom_point()+
   scale_colour_viridis_d(option="C")+
  labs(x=" FeOx (mg/g soil)", y="[MAOM N] (mg/g soil)", colour="Site")+
  theme_cowplot()+
  theme(legend.title=element_text(hjust=0.5, size=16), legend.text=element_text(size=15), axis.text=element_text(size=15), axis.title=element_text(size=19))

ggarrange(concC.fe.hf,concN.fe.hf, nrow=1, ncol=2, common.legend=T, legend="right")

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

ggsave("hf.conc.CN.fe.pdf",path="figures/", width= 20, height= 10, units="cm")

#what are the slopes of these lines?
fe_conchf_slopes_C=coef(lmList(mg.C.per.g.soil_HF~FeOx|Site , data = final_data))[2]
fe.mod1.slopes= as.data.frame(forests)%>% 
  arrange(forests) %>% 
  cbind(fe_conchf_slopes_C$FeOx) %>% 
  summarize(mean_slope=mean(fe_conchf_slopes_C$FeOx))

fe_conchf_slopes_N=coef(lmList(mg.N.per.g.soil_HF~FeOx|Site , data = final_data))[2]
fe.mod2.slopes= as.data.frame(forests)%>% 
  arrange(forests) %>% 
  cbind(fe_conchf_slopes_N$FeOx) %>% 
  summarize(mean_slope=mean(fe_conchf_slopes_N$FeOx))

##         1         2 
## -26.40569 -26.45323

	d13C					d14C
Predictors	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p
(Intercept)	-26.25	0.24	-0.00	0.14	<0.001	38.71	51.93	-0.01	0.40	0.461
FeOx	0.07	0.02	0.64	0.20	0.002	1.89	1.39	0.19	0.14	0.185
CDI	-17.29	7.45	-0.45	0.19	0.026	-1579.98	1489.10	-0.44	0.42	0.296
E	-0.03	0.21	-0.02	0.15	0.904	-12.25	13.00	-0.10	0.10	0.353
Random Effects
σ2	0.16					564.42
τ00	0.00 Site					1567.09 Site
ICC						0.74
N	6 Site					6 Site
Observations	40					40
Marginal R2 / Conditional R2	0.230 / NA					0.089 / 0.759

	d14C					d14C					d13C					d13C
Predictors	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p
(Intercept)	-12.86	14.42	0.00	0.16	0.378	-56.20	74.97	0.00	0.16	0.458	-26.42	0.16	-0.00	0.16	<0.001	-27.08	0.80	-0.00	0.16	<0.001
mg C per g soil HF	0.21	0.74	0.05	0.16	0.775						-0.00	0.01	-0.01	0.16	0.973
proportion C HF						58.01	92.10	0.10	0.16	0.533						0.81	0.99	0.13	0.16	0.419
Observations	40					40					40					40
R2 / R2 adjusted	0.002 / -0.024					0.010 / -0.016					0.000 / -0.026					0.017 / -0.009

	d14C					d14C					d14C					d14C					d14C					d14C
Predictors	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p
(Intercept)	-96.48	51.54	-0.00	0.42	0.135	-27.34	22.13	0.00	0.43	0.284	43.34	26.87	-0.00	0.33	0.182	48.99	65.67	-0.00	0.46	0.510	-118.57	64.15	0.00	0.35	0.138	-7.89	32.31	-0.00	0.43	0.823
FeOx	4.69	5.25	0.41	0.45	0.422	0.52	2.99	0.13	0.74	0.870	5.40	6.35	0.43	0.50	0.443	-10.42	11.74	-0.48	0.54	0.440	6.62	6.00	0.66	0.59	0.332	7.68	9.80	0.44	0.56	0.490
E	-22.78	49.88	-0.21	0.45	0.672	17.43	45.33	0.28	0.74	0.720	-17.01	26.68	-0.32	0.50	0.558	-26.97	46.72	-0.31	0.54	0.604	115.10	68.60	1.00	0.59	0.169	-36.78	29.43	-0.70	0.56	0.300
Observations	7					7					7					6					7					6
R2 / R2 adjusted	0.187 / -0.219					0.156 / -0.267					0.480 / 0.219					0.226 / -0.290					0.420 / 0.130					0.344 / -0.093

	d13C					d13C					d13C					d13C					d13C					d13C
Predictors	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p
(Intercept)	-26.31	0.40	-0.00	0.44	<0.001	-27.23	0.46	0.00	0.32	<0.001	-27.17	0.46	-0.00	0.16	<0.001	-26.87	0.67	0.00	0.51	<0.001	-26.48	0.50	0.00	0.29	<0.001	-27.25	0.50	-0.00	0.31	<0.001
FeOx	0.00	0.04	0.01	0.48	0.981	0.05	0.06	0.42	0.56	0.489	0.30	0.11	0.67	0.25	0.053	0.05	0.12	0.24	0.59	0.713	-0.03	0.05	-0.30	0.49	0.571	0.28	0.15	0.73	0.40	0.167
E	-0.23	0.38	-0.29	0.48	0.579	0.57	0.94	0.34	0.56	0.577	-0.63	0.46	-0.34	0.25	0.239	-0.01	0.48	-0.02	0.59	0.978	0.57	0.53	0.53	0.49	0.341	0.15	0.46	0.14	0.40	0.759
Observations	7					7					7					6					7					6
R2 / R2 adjusted	0.084 / -0.374					0.521 / 0.281					0.877 / 0.815					0.061 / -0.565					0.611 / 0.416					0.664 / 0.440

ecm_order=received_samples_plot %>% 
  select(plotID, order) %>% 
  distinct()

## Adding missing grouping variables: `Site`

pctAngio=ggplot(plotPctAngio %>%
                  rename("Angiosperm"="fractionAngiosperm") %>% 
                  mutate(Gymnosperm=1-Angiosperm,
                         plotID=as.character(plotID))%>% 
                           pivot_longer(c(Angiosperm, Gymnosperm), names_to="gymAng", values_to="dominance") %>% 
                  left_join(ecm_order) , aes(x=order, y=dominance, fill=gymAng))+
  geom_bar(position="stack", stat="identity", colour="black")+
  scale_fill_manual( values=c("gray90", "gray20"))+
  xlab("Study Plot")+
  ylab("Proportion of Basal Area")+
  facet_wrap(~Site, scales="free_x")+
  theme_cowplot()+
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(),strip.background = element_rect( fill="white"),legend.title=element_blank())

## Joining, by = "plotID"

ggsave("pctAngioGymno.jpg", path="figures/", width= 20, height= 16, units="cm")

BAbyMYC=ggplot(final_data %>% 
                 left_join(ecm_order), aes(x=order, y=plotBA))+
  geom_bar( stat="identity", colour="black")+
  xlab("Study plots in order of increasing ECM dominance")+
  ylab("Basal Area")+
  facet_wrap(~Site, scales="free_x")+
  theme_cowplot()+
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank(),strip.background = element_rect( fill="white"),legend.title=element_blank())

## Joining, by = c("plotID", "Site")

ggsave("BAmyc.jpg", path="figures/", width= 20, height= 16, units="cm")


ggplot(final_data, aes(x=fractionAngiosperm, y=proportion.C_HF, colour=Site))+
   geom_smooth(method="lm",size=2, se=F)+
   geom_point(size=3, alpha=0.5)+
  labs(x=" Angiosperm dominance ", y="MAOM C (proportion)", colour="Site")+
  theme_cowplot()+
  theme(legend.title=element_text(hjust=0.5, size=16), legend.text=element_text(size=15), axis.text=element_text(size=15), axis.title=element_text(size=19))

## `geom_smooth()` using formula 'y ~ x'

ggsave("maom.ang.prop.jpg", path="figures/", width= 19, height= 16, units="cm")

## `geom_smooth()` using formula 'y ~ x'

ggplot(final_data, aes(x=fractionAngiosperm, y=mg.C.per.g.soil_HF, colour=Site))+
   geom_smooth(method="lm",size=2, se=F)+
   geom_point(size=3, alpha=0.5)+
  labs(x=" Angiosperm dominance ", y="[MAOM C] (mg/g soil)", colour="Site")+
  theme_cowplot()+
  theme(legend.title=element_text(hjust=0.5, size=16), legend.text=element_text(size=15), axis.text=element_text(size=15), axis.title=element_text(size=19))

## `geom_smooth()` using formula 'y ~ x'

ggsave("maom.ang.conc.jpg", path="figures/", width= 19, height= 16, units="cm")

## `geom_smooth()` using formula 'y ~ x'

	mg.C.per.g.soil_HF					mg.N.per.g.soil_HF
Predictors	Estimates	std. Error	std. Beta	standardized std. Error	p	Estimates	std. Error	std. Beta	standardized std. Error	p
(Intercept)	15.364	3.487	0.017	0.285	<0.001	1.245	0.185	-0.002	0.204	<0.001
plotBA	0.539	1.329	0.055	0.135	0.687	-0.042	0.087	-0.074	0.154	0.632
Random Effects
σ2	51.014					0.251
τ00	40.908 Site					0.044 Site
ICC	0.445					0.150
N	7 Site					7 Site
Observations	47					47
Marginal R2 / Conditional R2	0.003 / 0.447					0.005 / 0.155

#Relationships between site-level climate decomposition index (CDI) and a) soil oxalate-extractable iron content, and b) total tree basal area in our study plots.

basal_area_cdi=ggplot(final_data, aes(x=CDI, y=plotBA))+
   geom_point()+
  stat_summary(fun= mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.0009, position="dodge")+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.0005, position = position_dodge(width = 0.0009))+
                theme_cowplot()+
  xlab("Climate Decomposition Index")+
  ylab("Total Basal Area"
  )


feox_cdi=ggplot(final_data, aes(x=CDI, y=FeOx))+
   geom_point()+
  stat_summary(fun= mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.0009, position="dodge")+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.0005, position = position_dodge(width = 0.0009))+
    geom_smooth(method="lm")+
                theme_cowplot()+
  xlab("Climate Decomposition Index")+
  ylab("Soil Oxalate-Extractable Iron"
  )

ggarrange(basal_area_cdi, feox_cdi, ncol=2, nrow=1, labels= c("a","b"), legend=F)

## `geom_smooth()` using formula 'y ~ x'

ggsave("FigS2.pdf", path="figures/", width= 24, height= 12, units="cm")

summary(lm(plotBA~CDI,data=final_data )) #

## 
## Call:
## lm(formula = plotBA ~ CDI, data = final_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7455 -0.7491 -0.1379  0.6662  2.1157 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.7240     0.4667   3.694 0.000596 ***
## CDI           0.5506    12.8459   0.043 0.966001    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9489 on 45 degrees of freedom
## Multiple R-squared:  4.082e-05,  Adjusted R-squared:  -0.02218 
## F-statistic: 0.001837 on 1 and 45 DF,  p-value: 0.966

summary(lm(FeOx~CDI,data=final_data )) #yes, more FeOx in warmer sites

## 
## Call:
## lm(formula = FeOx ~ CDI, data = final_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.5624 -2.8296  0.2597  1.6610 11.9572 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -1.093      1.821  -0.600 0.551306    
## CDI          184.539     50.119   3.682 0.000617 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.702 on 45 degrees of freedom
## Multiple R-squared:  0.2315, Adjusted R-squared:  0.2144 
## F-statistic: 13.56 on 1 and 45 DF,  p-value: 0.0006175

summary(lm(mg.C.per.g.soil_HF~FeOx,data=final_data %>% filter(Site=="LENO") )) #yes within this site Fe Ox is a big driver of concentrations

## 
## Call:
## lm(formula = mg.C.per.g.soil_HF ~ FeOx, data = final_data %>% 
##     filter(Site == "LENO"))
## 
## Residuals:
##       1       2       3       4       5       6       7 
## -4.0951  4.7057 -1.6980  2.2892 -1.6658  0.1938  0.2702 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   3.5803     2.5703   1.393  0.22239   
## FeOx          0.9075     0.2157   4.207  0.00843 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.16 on 5 degrees of freedom
## Multiple R-squared:  0.7798, Adjusted R-squared:  0.7357 
## F-statistic:  17.7 on 1 and 5 DF,  p-value: 0.008429

summary(lm(mg.C.per.g.soil_HF~FeOx,data=final_data %>% filter(Site=="DELA") ))  #yes within this site Fe Ox is a big driver of concentrations

## 
## Call:
## lm(formula = mg.C.per.g.soil_HF ~ FeOx, data = final_data %>% 
##     filter(Site == "DELA"))
## 
## Residuals:
##       1       2       3       4       5       6       7 
## -1.8847 -0.2968  3.2730  3.2409 -4.2195  1.1282 -1.2411 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   2.6878     4.0998   0.656   0.5410  
## FeOx          1.1080     0.4278   2.590   0.0488 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.016 on 5 degrees of freedom
## Multiple R-squared:  0.5729, Adjusted R-squared:  0.4875 
## F-statistic: 6.707 on 1 and 5 DF,  p-value: 0.04884

summary(lm(FeOx~E,data=final_data )) #no overall trend with FeOx and ECM dominance across sites

## 
## Call:
## lm(formula = FeOx ~ E, data = final_data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.211 -2.730 -0.970  2.103 13.966 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.182      1.125   3.716 0.000557 ***
## E              2.301      1.933   1.190 0.240214    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.158 on 45 degrees of freedom
## Multiple R-squared:  0.03052,    Adjusted R-squared:  0.008974 
## F-statistic: 1.417 on 1 and 45 DF,  p-value: 0.2402

summary(lm(FeOx~E,data=final_data %>% filter(Site=="LENO") )) # yes within LENO positive relationship between soil FeOx and ECM dominance

## 
## Call:
## lm(formula = FeOx ~ E, data = final_data %>% filter(Site == "LENO"))
## 
## Residuals:
##       1       2       3       4       5       6       7 
##  6.4818 -0.3810 -0.9449 -2.4064 -1.4924  3.2959 -4.5530 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)    4.001      2.782   1.438   0.2099  
## E             11.881      4.206   2.825   0.0369 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.066 on 5 degrees of freedom
## Multiple R-squared:  0.6148, Adjusted R-squared:  0.5378 
## F-statistic: 7.981 on 1 and 5 DF,  p-value: 0.03689

summary(lm(FeOx~E,data=final_data %>% filter(Site=="DELA") )) # not at all within DELA 

## 
## Call:
## lm(formula = FeOx ~ E, data = final_data %>% filter(Site == "DELA"))
## 
## Residuals:
##       1       2       3       4       5       6       7 
## -5.2936 -0.5310 -0.5159 -0.7887  1.2474  2.2708  3.6111 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)    8.768      1.978   4.433  0.00681 **
## E              1.162      4.219   0.275  0.79401   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.129 on 5 degrees of freedom
## Multiple R-squared:  0.01495,    Adjusted R-squared:  -0.1821 
## F-statistic: 0.07586 on 1 and 5 DF,  p-value: 0.794

#Site means of MAOM C and N concentrations and proportions, ordered by climate decomposition index

MAOM_C_conc_cdi=ggplot(final_data, aes(x=CDI, y=mg.C.per.g.soil_HF))+
   geom_point()+
  stat_summary(fun= mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.0009, position="dodge")+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.0005, position = position_dodge(width = 0.0009))+
  geom_smooth(method="lm")+
                theme_cowplot()+
  xlab("Climate Decomposition Index")+
  ylab("[MAOM C] (mg/g soil)"
  )

MAOM_C_prop_cdi=ggplot(final_data, aes(x=CDI, y=proportion.C_HF))+
   geom_point()+
  stat_summary(fun= mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.0009, position="dodge")+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.0005, position = position_dodge(width = 0.0009))+
  theme_cowplot()+
  xlab("Climate Decomposition Index")+
  ylab(expression("proportion C"[MAOM])
  )

MAOM_N_conc_cdi=ggplot(final_data, aes(x=CDI, y=mg.N.per.g.soil_HF))+
   geom_point()+
  stat_summary(fun= mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.0009, position="dodge")+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.0005, position = position_dodge(width = 0.0009))+
  theme_cowplot()+
  xlab("Climate Decomposition Index")+
  ylab("[MAOM N] (mg/g soil)"
  )

MAOM_N_prop_cdi=ggplot(final_data, aes(x=CDI, y=proportion.N_HF))+
   geom_point()+
  stat_summary(fun= mean, fun.min=mean, fun.max=mean, geom="crossbar", width=0.0009, position="dodge")+
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.0005, position = position_dodge(width = 0.0009))+
  geom_smooth(method="lm")+
                theme_cowplot()+
  xlab("Climate Decomposition Index")+
  ylab(expression("proportion N"[MAOM])
  )

ggarrange(MAOM_C_conc_cdi, MAOM_C_prop_cdi, MAOM_N_conc_cdi, MAOM_N_prop_cdi, ncol=2, nrow=2, labels= c("a","b", "c", "d"), legend=F)

## `geom_smooth()` using formula 'y ~ x'
## `geom_smooth()` using formula 'y ~ x'

ggsave("FigXXX.jpg", path="figures/", width= 24, height= 22, units="cm")

#Is the amount of oxalate-extractable iron in soil related to mycorrhizal type?

ggplot(final_data, aes(x=E, y=FeOx))+
  geom_point(color="blue")+
  geom_smooth(method="lm", color="blue")+
  facet_wrap(~Site)+
  theme_cowplot()

## `geom_smooth()` using formula 'y ~ x'

ggsave("ECM_vs_FeOx.jpg", path="figures/", width= 30, height= 24, units="cm")

## `geom_smooth()` using formula 'y ~ x'

ggplot(final_data, aes(x=CDI, y=FeOx))+
  geom_point(color="blue")+
  geom_smooth(method="lm", color="blue")+
  theme_cowplot()

## `geom_smooth()` using formula 'y ~ x'

ggsave("CDI_vs_FeOx.pdf", path="figures/", width= 30, height= 24, units="cm")

## `geom_smooth()` using formula 'y ~ x'

	C.N_FLF						C.N_OLF						C.N_HF
Predictors	Estimates	std. Error	std. Beta	standardized std. Error	p	df	Estimates	std. Error	std. Beta	standardized std. Error	p	df	Estimates	std. Error	std. Beta	standardized std. Error	p	df
(Intercept)	23.724	3.861	0.000	0.174	<0.001	41.000	22.582	5.534	-0.001	0.199	<0.001	41.000	23.195	2.405	-0.001	0.160	<0.001	41.000
E	7.898	2.508	0.413	0.131	0.003	41.000	12.146	3.133	0.494	0.127	<0.001	41.000	2.322	1.110	0.171	0.082	0.043	41.000
CDI	66.397	108.571	0.119	0.195	0.544	41.000	141.022	156.075	0.197	0.218	0.372	41.000	-299.874	67.826	-0.757	0.171	<0.001	41.000
FeOx	-0.641	0.245	-0.441	0.168	0.012	41.000	-0.428	0.319	-0.229	0.171	0.187	41.000	-0.013	0.118	-0.013	0.114	0.912	41.000
Random Effects
σ2	25.801						39.696						4.909
τ00	3.948 Site						10.891 Site						2.607 Site
ICC	0.133						0.215						0.347
N	7 Site						7 Site						7 Site
Observations	47						47						47
Marginal R2 / Conditional R2	0.241 / 0.341						0.220 / 0.388						0.611 / 0.746

	C.N_HF		C.N_HF		C.N_HF		C.N_HF		C.N_HF		C.N_HF		C.N_HF
Predictors	Estimates	p	Estimates	p	Estimates	p	Estimates	p	Estimates	p	Estimates	p	Estimates	p
(Intercept)	20.25	0.001	16.55	0.002	13.96	0.003	10.98	<0.001	10.58	0.006	10.00	<0.001	9.77	<0.001
E	-2.03	0.701	3.25	0.347	4.18	0.299	3.25	0.221	4.71	0.331	-0.46	0.540	-0.76	0.183
Observations	7		6		6		7		7		7		7
R2 / R2 adjusted	0.032 / -0.162		0.221 / 0.026		0.262 / 0.078		0.281 / 0.138		0.188 / 0.025		0.080 / -0.104		0.323 / 0.188

`